{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "HelloGPU"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "void helloCPU()\n",
    "{\n",
    "  printf(\"Hello from the CPU.\\n\");\n",
    "}\n",
    "\n",
    "/*\n",
    " * Refactor the `helloGPU` definition to be a kernel\n",
    " * that can be launched on the GPU. Update its message\n",
    " * to read \"Hello from the GPU!\"\n",
    " */\n",
    "\n",
    "__global__ void helloGPU()\n",
    "{\n",
    "  printf(\"Hello from the GPU.\\n\");\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  helloGPU()<<<1,1>>>;\n",
    "  helloCPU();\n",
    "\n",
    "  /*\n",
    "   * Refactor this call to `helloGPU` so that it launches\n",
    "   * as a kernel on the GPU.\n",
    "   */\n",
    "\n",
    "  helloGPU()<<<1,1>>>;\n",
    "\n",
    "  /*\n",
    "   * Add code below to synchronize on the completion of the\n",
    "   * `helloGPU` kernel completion before continuing the CPU\n",
    "   * thread.\n",
    "   */\n",
    "\n",
    "   cudaDeviceSynchronize();\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "\n",
    "/*\n",
    " * Refactor firstParallel so that it can run on the GPU.\n",
    " */\n",
    "\n",
    "__global__ void firstParallel()\n",
    "{\n",
    "  printf(\"This should be running in parallel.\\n\");\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  /*\n",
    "   * Refactor this call to firstParallel to execute in parallel\n",
    "   * on the GPU.\n",
    "   */\n",
    "\n",
    "  firstParallel<<<1,5>>>();\n",
    "  cudaDeviceSynchronize();\n",
    "\n",
    "  /*\n",
    "   * Some code is needed below so that the CPU will wait\n",
    "   * for the GPU kernels to complete before proceeding.\n",
    "   */\n",
    "\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "使用特点的块和线程索引"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "\n",
    "__global__ void printSuccessForCorrectExecutionConfiguration()\n",
    "{\n",
    "\n",
    "  if(threadIdx.x == 1023 && blockIdx.x == 255)\n",
    "  {\n",
    "    printf(\"Success!\\n\");\n",
    "  } else {\n",
    "    printf(\"Failure. Update the execution configuration as necessary.\\n\");\n",
    "  }\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  /*\n",
    "   * Update the execution configuration so that the kernel\n",
    "   * will print `\"Success!\"`.\n",
    "   */\n",
    "\n",
    "  printSuccessForCorrectExecutionConfiguration<<<256, 1024>>>();\n",
    "  cudaDeviceSynchronize();\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "加速单块for循环"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "\n",
    "/*\n",
    " * Refactor `loop` to be a CUDA Kernel. The new kernel should\n",
    " * only do the work of 1 iteration of the original loop.\n",
    " */\n",
    "\n",
    "__global__ void loop(int N)\n",
    "{\n",
    "  for (int i = threadIdx.x*2; i < (threadIdx.x+1)*2; ++i)\n",
    "  {\n",
    "    printf(\"This is iteration number %d\\n\", i);\n",
    "  }\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  /*\n",
    "   * When refactoring `loop` to launch as a kernel, be sure\n",
    "   * to use the execution configuration to control how many\n",
    "   * \"iterations\" to perform.\n",
    "   *\n",
    "   * For this exercise, only use 1 block of threads.\n",
    "   */\n",
    "\n",
    "  int N = 10;\n",
    "  loop<<<1,5>>>(N);\n",
    "  cudaDeviceSynchronize();\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "加速更多数字的for循环。与上例不同，这里需要启动更多线程块，每个线程块负责一定的任务"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "\n",
    "/*\n",
    " * Refactor `loop` to be a CUDA Kernel. The new kernel should\n",
    " * only do the work of 1 iteration of the original loop.\n",
    " */\n",
    "\n",
    "__global__ void loop(int N)\n",
    "{\n",
    "    printf(\"%d\\n\",threadIdx.x+blockIdx.x*blockDim.x);\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  /*\n",
    "   * When refactoring `loop` to launch as a kernel, be sure\n",
    "   * to use the execution configuration to control how many\n",
    "   * \"iterations\" to perform.\n",
    "   *\n",
    "   * For this exercise, be sure to use more than 1 block in\n",
    "   * the execution configuration.\n",
    "   */\n",
    "\n",
    "  int N = 10;\n",
    "  loop<<<10,10>>>(N);\n",
    "  cudaDeviceSynchronize();\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "内存在CPU和GPU上的分配"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "\n",
    "/*\n",
    " * Initialize array values on the host.\n",
    " */\n",
    "\n",
    "void init(int *a, int N)\n",
    "{\n",
    "  int i;\n",
    "  for (i = 0; i < N; ++i)\n",
    "  {\n",
    "    a[i] = i;\n",
    "  }\n",
    "}\n",
    "\n",
    "/*\n",
    " * Double elements in parallel on the GPU.\n",
    " */\n",
    "\n",
    "__global__\n",
    "void doubleElements(int *a, int N)\n",
    "{\n",
    "  int i;\n",
    "  i = blockIdx.x * blockDim.x + threadIdx.x;\n",
    "  if (i < N)\n",
    "  {\n",
    "    a[i] *= 2;\n",
    "  }\n",
    "}\n",
    "\n",
    "/*\n",
    " * Check all elements have been doubled on the host.\n",
    " */\n",
    "\n",
    "bool checkElementsAreDoubled(int *a, int N)\n",
    "{\n",
    "  int i;\n",
    "  for (i = 0; i < N; ++i)\n",
    "  {\n",
    "    if (a[i] != i*2) return false;\n",
    "  }\n",
    "  return true;\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  int N = 100;\n",
    "  int *a;\n",
    "\n",
    "  size_t size = N * sizeof(int);\n",
    "\n",
    "  /*\n",
    "   * Refactor this memory allocation to provide a pointer\n",
    "   * `a` that can be used on both the host and the device.\n",
    "   */\n",
    "\n",
    "  cudaMallocManaged(&a, size);\n",
    "  //a = (int *)malloc(size);\n",
    "\n",
    "  init(a, N);\n",
    "\n",
    "  size_t threads_per_block = 10;\n",
    "  size_t number_of_blocks = 10;\n",
    "\n",
    "  /*\n",
    "   * This launch will not work until the pointer `a` is also\n",
    "   * available to the device.\n",
    "   */\n",
    "\n",
    "  doubleElements<<<number_of_blocks, threads_per_block>>>(a, N);\n",
    "  cudaDeviceSynchronize();\n",
    "\n",
    "  bool areDoubled = checkElementsAreDoubled(a, N);\n",
    "  printf(\"All elements were doubled? %s\\n\", areDoubled ? \"TRUE\" : \"FALSE\");\n",
    "\n",
    "  /*\n",
    "   * Refactor to free memory that has been allocated to be\n",
    "   * accessed by both the host and the device.\n",
    "   */\n",
    "\n",
    "  cudaFree(a);\n",
    "  //free(a);\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "当线程数与任务量n不匹配时，一种思路是创建足够多的线程（32的倍数），对每个线程计算其各自的id = threadIdx.x+blockIdx.x*blockDim.x，若id < n，则执行任务，否则闲置"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "\n",
    "/*\n",
    " * Currently, `initializeElementsTo`, if executed in a thread whose\n",
    " * `i` is calculated to be greater than `N`, will try to access a value\n",
    " * outside the range of `a`.\n",
    " *\n",
    " * Refactor the kernel definition to prevent out of range accesses.\n",
    " */\n",
    "\n",
    "__global__ void initializeElementsTo(int initialValue, int *a, int N)\n",
    "{\n",
    "  int i = threadIdx.x + blockIdx.x * blockDim.x;\n",
    "  if(i<N)   a[i] = initialValue;\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  /*\n",
    "   * Do not modify `N`.\n",
    "   */\n",
    "\n",
    "  int N = 1000;\n",
    "\n",
    "  int *a;\n",
    "  size_t size = N * sizeof(int);\n",
    "\n",
    "  cudaMallocManaged(&a, size);\n",
    "\n",
    "  /*\n",
    "   * Assume we have reason to want the number of threads\n",
    "   * fixed at `256`: do not modify `threads_per_block`.\n",
    "   */\n",
    "\n",
    "  size_t threads_per_block = 256;\n",
    "\n",
    "  /*\n",
    "   * Assign a value to `number_of_blocks` that will\n",
    "   * allow for a working execution configuration given\n",
    "   * the fixed values for `N` and `threads_per_block`.\n",
    "   */\n",
    "\n",
    "  size_t number_of_blocks = 4;\n",
    "\n",
    "  int initialValue = 6;\n",
    "\n",
    "  initializeElementsTo<<<number_of_blocks, threads_per_block>>>(initialValue, a, N);\n",
    "  cudaDeviceSynchronize();\n",
    "\n",
    "  /*\n",
    "   * Check to make sure all values in `a`, were initialized.\n",
    "   */\n",
    "\n",
    "  for (int i = 0; i < N; ++i)\n",
    "  {\n",
    "    if(a[i] != initialValue)\n",
    "    {\n",
    "      printf(\"FAILURE: target value: %d\\t a[%d]: %d\\n\", initialValue, i, a[i]);\n",
    "      cudaFree(a);\n",
    "      exit(1);\n",
    "    }\n",
    "  }\n",
    "  printf(\"SUCCESS!\\n\");\n",
    "\n",
    "  cudaFree(a);\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "另一种思路，是创建线程数少于任务量的grid，通过grid-stride循环获取对应的任务，每个线程将计算其在网格内的唯一索引，对数组中该索引处的元素执行其操作，然后将索引增加stride，直到它超出任务的范围。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "\n",
    "void init(int *a, int N)\n",
    "{\n",
    "  int i;\n",
    "  for (i = 0; i < N; ++i)\n",
    "  {\n",
    "    a[i] = i;\n",
    "  }\n",
    "}\n",
    "\n",
    "/*\n",
    " * In the current application, `N` is larger than the grid.\n",
    " * Refactor this kernel to use a grid-stride loop in order that\n",
    " * each parallel thread work on more than one element of the array.\n",
    " */\n",
    "\n",
    "__global__\n",
    "void doubleElements(int *a, int N)\n",
    "{\n",
    "  int i, grid;\n",
    "  i = blockIdx.x * blockDim.x + threadIdx.x;\n",
    "  grid = blockDim.x * gridDim.x;\n",
    "  while (i < N)\n",
    "  {\n",
    "    a[i] *= 2;\n",
    "    i += grid;\n",
    "  }\n",
    "}\n",
    "\n",
    "bool checkElementsAreDoubled(int *a, int N)\n",
    "{\n",
    "  int i;\n",
    "  for (i = 0; i < N; ++i)\n",
    "  {\n",
    "    if (a[i] != i*2) return false;\n",
    "  }\n",
    "  return true;\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  /*\n",
    "   * `N` is greater than the size of the grid (see below).\n",
    "   */\n",
    "\n",
    "  int N = 10000;\n",
    "  int *a;\n",
    "\n",
    "  size_t size = N * sizeof(int);\n",
    "  cudaMallocManaged(&a, size);\n",
    "\n",
    "  init(a, N);\n",
    "\n",
    "  /*\n",
    "   * The size of this grid is 256*32 = 8192.\n",
    "   */\n",
    "\n",
    "  size_t threads_per_block = 256;\n",
    "  size_t number_of_blocks = 32;\n",
    "\n",
    "  doubleElements<<<number_of_blocks, threads_per_block>>>(a, N);\n",
    "  cudaDeviceSynchronize();\n",
    "\n",
    "  bool areDoubled = checkElementsAreDoubled(a, N);\n",
    "  printf(\"All elements were doubled? %s\\n\", areDoubled ? \"TRUE\" : \"FALSE\");\n",
    "\n",
    "  cudaFree(a);\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "错误类型"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "#include <assert.h>\n",
    "\n",
    "inline cudaError_t checkCuda(cudaError_t result)\n",
    "{\n",
    "  if (result != cudaSuccess) {\n",
    "    fprintf(stderr, \"CUDA Runtime Error: %s\\n\", cudaGetErrorString(result));\n",
    "    assert(result == cudaSuccess);\n",
    "  }\n",
    "  return result;\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "\n",
    "/*\n",
    " * The macro can be wrapped around any function returning\n",
    " * a value of type `cudaError_t`.\n",
    " */\n",
    "\n",
    "  checkCuda( cudaDeviceSynchronize() )\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Final Exercise: Accelerate Vector Addition Application"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "1. Augment the addVectorsInto definition so that it is a CUDA kernel.\n",
    "2. Choose and utilize a working execution configuration so that addVectorsInto launches as a CUDA kernel.\n",
    "3. Update memory allocations, and memory freeing to reflect that the 3 vectors a, b, and result need to be accessed by host and device code.\n",
    "4. Refactor the body of addVectorsInto: it will be launched inside of a single thread, and only needs to do one thread's worth of work on the input vectors. Be certain the thread will never try to access elements outside the range of the input vectors, and take care to note whether or not the thread needs to do work on more than one element of the input vectors.\n",
    "5. Add error handling in locations where CUDA code might otherwise silently fail."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "#include <assert.h>\n",
    "\n",
    "inline cudaError_t checkCuda(cudaError_t result)\n",
    "{\n",
    "  if (result != cudaSuccess) {\n",
    "    fprintf(stderr, \"CUDA Runtime Error: %s\\n\", cudaGetErrorString(result));\n",
    "    assert(result == cudaSuccess);\n",
    "  }\n",
    "  return result;\n",
    "}\n",
    "\n",
    "void initWith(float num, float *a, int N)\n",
    "{\n",
    "  for(int i = 0; i < N; ++i)\n",
    "  {\n",
    "    a[i] = num;\n",
    "  }\n",
    "}\n",
    "\n",
    "__global__ void addVectorsInto(float *result, float *a, float *b, int N)\n",
    "{\n",
    "  int id = threadIdx.x +blockIdx.x * blockDim.x;\n",
    "  int stride = blockDim.x * gridDim.x;\n",
    "  while(id < N)\n",
    "  {\n",
    "    result[id] = a[id] + b[id];\n",
    "    id += stride;\n",
    "  }\n",
    "}\n",
    "\n",
    "void checkElementsAre(float target, float *array, int N)\n",
    "{\n",
    "  for(int i = 0; i < N; i++)\n",
    "  {\n",
    "    if(array[i] != target)\n",
    "    {\n",
    "      printf(\"FAIL: array[%d] - %0.0f does not equal %0.0f\\n\", i, array[i], target);\n",
    "      exit(1);\n",
    "    }\n",
    "  }\n",
    "  printf(\"SUCCESS! All values added correctly.\\n\");\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  const int N = 2<<20;\n",
    "  size_t size = N * sizeof(float);\n",
    "\n",
    "  size_t threads_per_block = 1024;\n",
    "  size_t blocks_per_grid = (N + threads_per_block - 1) / threads_per_block;\n",
    "\n",
    "  float *a;\n",
    "  float *b;\n",
    "  float *c;\n",
    "\n",
    "  checkCuda(cudaMallocManaged(&a, size));\n",
    "  checkCuda(cudaMallocManaged(&b, size));\n",
    "  checkCuda(cudaMallocManaged(&c, size));\n",
    "\n",
    "  initWith(3, a, N);\n",
    "  initWith(4, b, N);\n",
    "  initWith(0, c, N);\n",
    "\n",
    "  addVectorsInto<<<blocks_per_grid,threads_per_block>>>(c, a, b, N);\n",
    "  checkCuda(cudaDeviceSynchronize());\n",
    "  checkCuda(cudaGetLastError());\n",
    "\n",
    "  checkElementsAre(7, c, N);\n",
    "\n",
    "  checkCuda(cudaFree(a));\n",
    "  checkCuda(cudaFree(b));\n",
    "  checkCuda(cudaFree(c));\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "2或3维的grid和block，使用dim3 name(x_size,y_size,z_size);\n",
    "\n",
    "相应block的grid的规模可以分别通过 dim3.x 和 dim3.y 获取"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accelerate 2D Matrix Multiply Application\n",
    "1. You will need to create an execution configuration whose arguments are both values with the and dimensions set to greater than .dim3xy1\n",
    "2. Inside the body of the kernel, you will need to establish the running thread's unique index within the grid per usual, but you should establish two indices for the thread: one for the x axis of the grid, and one for the y axis of the grid."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "\n",
    "#define N  64\n",
    "\n",
    "__global__ void matrixMulGPU( int * a, int * b, int * c )\n",
    "{\n",
    "  /*\n",
    "   * Build out this kernel.\n",
    "   */\n",
    "   int value = 0;\n",
    "   int row = threadIdx.x + blockDim.x * blockIdx.x;\n",
    "   int col = threadIdx.y + blockDim.y * blockIdx.y;\n",
    "\n",
    "   if(row < N && col < N)\n",
    "   {\n",
    "    for(int k = 0; k < N; ++k)\n",
    "    {\n",
    "      value += a[row * N + k] * b[k * N + col];\n",
    "    }\n",
    "      c[row * N + col] = value;\n",
    "   }\n",
    "}\n",
    "\n",
    "/*\n",
    " * This CPU function already works, and will run to create a solution matrix\n",
    " * against which to verify your work building out the matrixMulGPU kernel.\n",
    " */\n",
    "\n",
    "void matrixMulCPU( int * a, int * b, int * c )\n",
    "{\n",
    "  int val = 0;\n",
    "\n",
    "  for( int row = 0; row < N; ++row )\n",
    "    for( int col = 0; col < N; ++col )\n",
    "    {\n",
    "      val = 0;\n",
    "      for ( int k = 0; k < N; ++k )\n",
    "        val += a[row * N + k] * b[k * N + col];\n",
    "      c[row * N + col] = val;\n",
    "    }\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  int *a, *b, *c_cpu, *c_gpu; // Allocate a solution matrix for both the CPU and the GPU operations\n",
    "\n",
    "  int size = N * N * sizeof (int); // Number of bytes of an N x N matrix\n",
    "\n",
    "  // Allocate memory\n",
    "  cudaMallocManaged (&a, size);\n",
    "  cudaMallocManaged (&b, size);\n",
    "  cudaMallocManaged (&c_cpu, size);\n",
    "  cudaMallocManaged (&c_gpu, size);\n",
    "\n",
    "  // Initialize memory; create 2D matrices\n",
    "  for( int row = 0; row < N; ++row )\n",
    "    for( int col = 0; col < N; ++col )\n",
    "    {\n",
    "      a[row*N + col] = row;\n",
    "      b[row*N + col] = col+2;\n",
    "      c_cpu[row*N + col] = 0;\n",
    "      c_gpu[row*N + col] = 0;\n",
    "    }\n",
    "\n",
    "  /*\n",
    "   * Assign `threads_per_block` and `number_of_blocks` 2D values\n",
    "   * that can be used in matrixMulGPU above.\n",
    "   */\n",
    "\n",
    "  dim3 threads_per_block(16,16,1);\n",
    "  dim3 number_of_blocks(N/threads_per_block.x+1,N/threads_per_block.y+1,1);\n",
    "\n",
    "  matrixMulGPU <<< number_of_blocks, threads_per_block >>> ( a, b, c_gpu );\n",
    "\n",
    "  cudaDeviceSynchronize();\n",
    "\n",
    "  // Call the CPU version to check our work\n",
    "  matrixMulCPU( a, b, c_cpu );\n",
    "\n",
    "  // Compare the two answers to make sure they are equal\n",
    "  bool error = false;\n",
    "  for( int row = 0; row < N && !error; ++row )\n",
    "    for( int col = 0; col < N && !error; ++col )\n",
    "      if (c_cpu[row * N + col] != c_gpu[row * N + col])\n",
    "      {\n",
    "        printf(\"FOUND ERROR at c[%d][%d]\\n\", row, col);\n",
    "        error = true;\n",
    "        break;\n",
    "      }\n",
    "  if (!error)\n",
    "    printf(\"Success!\\n\");\n",
    "\n",
    "  // Free all our allocated memory\n",
    "  cudaFree(a); cudaFree(b);\n",
    "  cudaFree( c_cpu ); cudaFree( c_gpu );\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Accelerate A Thermal Conductivity Application\n",
    "\n",
    "Convert the step_kernel_mod function inside 01-heat-conduction.cu to execute on the GPU, and modify the main function to properly allocate data for use on CPU and GPU. The step_kernel_ref function executes on the CPU and is used for error checking. Because this code involves floating point calculations, different processors, or even simply reordering operations on the same processor, can result in slightly different results. For this reason the error checking code uses an error threshold, instead of looking for an exact match.\n",
    "\n",
    "浮点计算结果随设备环境不同，这里的错误判断就依据一定的阈值"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#include <stdio.h>\n",
    "#include <math.h>\n",
    "\n",
    "// Simple define to index into a 1D array from 2D space\n",
    "#define I2D(num, c, r) ((r)*(num)+(c))\n",
    "\n",
    "/*\n",
    " * `step_kernel_mod` is currently a direct copy of the CPU reference solution\n",
    " * `step_kernel_ref` below. Accelerate it to run as a CUDA kernel.\n",
    " */\n",
    "\n",
    "__global__ \n",
    "void step_kernel_mod(int ni, int nj, float fact, float* temp_in, float* temp_out)\n",
    "{\n",
    "  int i00, im10, ip10, i0m1, i0p1;\n",
    "  float d2tdx2, d2tdy2;\n",
    "\n",
    "  int j = threadIdx.x + blockDim.x * blockIdx.x;\n",
    "  int i = threadIdx.y + blockDim.y * blockIdx.y;\n",
    "\n",
    "  if((j >0 && j < nj-1)&&(i > 0 && i < ni - 1))\n",
    "  {\n",
    "      i00 = I2D(ni, i, j);\n",
    "      im10 = I2D(ni, i-1, j);\n",
    "      ip10 = I2D(ni, i+1, j);\n",
    "      i0m1 = I2D(ni, i, j-1);\n",
    "      i0p1 = I2D(ni, i, j+1);\n",
    "\n",
    "      d2tdx2 = temp_in[im10]-2*temp_in[i00]+temp_in[ip10];\n",
    "      d2tdy2 = temp_in[i0m1]-2*temp_in[i00]+temp_in[i0p1];\n",
    "\n",
    "      temp_out[i00] = temp_in[i00]+fact*(d2tdx2 + d2tdy2);\n",
    "  }\n",
    "}\n",
    "\n",
    "void step_kernel_ref(int ni, int nj, float fact, float* temp_in, float* temp_out)\n",
    "{\n",
    "  int i00, im10, ip10, i0m1, i0p1;\n",
    "  float d2tdx2, d2tdy2;\n",
    "\n",
    "\n",
    "  // loop over all points in domain (except boundary)\n",
    "  for ( int j=1; j < nj-1; j++ ) {\n",
    "    for ( int i=1; i < ni-1; i++ ) {\n",
    "      // find indices into linear memory\n",
    "      // for central point and neighbours\n",
    "      i00 = I2D(ni, i, j);\n",
    "      im10 = I2D(ni, i-1, j);\n",
    "      ip10 = I2D(ni, i+1, j);\n",
    "      i0m1 = I2D(ni, i, j-1);\n",
    "      i0p1 = I2D(ni, i, j+1);\n",
    "\n",
    "      // evaluate derivatives\n",
    "      d2tdx2 = temp_in[im10]-2*temp_in[i00]+temp_in[ip10];\n",
    "      d2tdy2 = temp_in[i0m1]-2*temp_in[i00]+temp_in[i0p1];\n",
    "\n",
    "      // update temperatures\n",
    "      temp_out[i00] = temp_in[i00]+fact*(d2tdx2 + d2tdy2);\n",
    "    }\n",
    "  }\n",
    "}\n",
    "\n",
    "int main()\n",
    "{\n",
    "  int istep;\n",
    "  int nstep = 200; // number of time steps\n",
    "\n",
    "  // Specify our 2D dimensions\n",
    "  const int ni = 200;\n",
    "  const int nj = 100;\n",
    "  float tfac = 8.418e-5; // thermal diffusivity of silver\n",
    "\n",
    "  float *temp1_ref, *temp2_ref, *temp1, *temp2, *temp_tmp;\n",
    "\n",
    "  const int size = ni * nj * sizeof(float);\n",
    "\n",
    "  temp1_ref = (float*)malloc(size);\n",
    "  temp2_ref = (float*)malloc(size);\n",
    "  cudaMallocManaged(&temp1, size);\n",
    "  cudaMallocManaged(&temp2, size);\n",
    "\n",
    "  // Initialize with random data\n",
    "  for( int i = 0; i < ni*nj; ++i) {\n",
    "    temp1_ref[i] = temp2_ref[i] = temp1[i] = temp2[i] = (float)rand()/(float)(RAND_MAX/100.0f);\n",
    "  }\n",
    "\n",
    "  // Execute the CPU-only reference version\n",
    "  for (istep=0; istep < nstep; istep++) {\n",
    "    step_kernel_ref(ni, nj, tfac, temp1_ref, temp2_ref);\n",
    "\n",
    "    // swap the temperature pointers\n",
    "    temp_tmp = temp1_ref;\n",
    "    temp1_ref = temp2_ref;\n",
    "    temp2_ref= temp_tmp;\n",
    "  }\n",
    "\n",
    "  // Execute the modified version using same data\n",
    "  dim3 threads_per_block(16,16,1);\n",
    "  dim3 blocks_per_grid(nj/threads_per_block.x+1,ni/threads_per_block.y+1,1);\n",
    "  for (istep=0; istep < nstep; istep++) {\n",
    "    step_kernel_mod<<<threads_per_block,blocks_per_grid>>>(ni, nj, tfac, temp1, temp2);\n",
    "    cudaDeviceSynchronize();\n",
    "\n",
    "    // swap the temperature pointers\n",
    "    temp_tmp = temp1;\n",
    "    temp1 = temp2;\n",
    "    temp2= temp_tmp;\n",
    "  }\n",
    "\n",
    "  float maxError = 0;\n",
    "  // Output should always be stored in the temp1 and temp1_ref at this point\n",
    "  for( int i = 0; i < ni*nj; ++i ) {\n",
    "    if (abs(temp1[i]-temp1_ref[i]) > maxError) { maxError = abs(temp1[i]-temp1_ref[i]); }\n",
    "  }\n",
    "\n",
    "  // Check and see if our maxError is greater than an error bound\n",
    "  if (maxError > 0.0005f)\n",
    "    printf(\"Problem! The Max Error of %.5f is NOT within acceptable bounds.\\n\", maxError);\n",
    "  else\n",
    "    printf(\"The Max Error of %.5f is within acceptable bounds.\\n\", maxError);\n",
    "\n",
    "free(temp1_ref);\n",
    "free(temp2_ref);\n",
    "cudaFree(temp1);\n",
    "cudaFree(temp2);\n",
    "  return 0;\n",
    "}\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3.10.2 64-bit",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "name": "python",
   "version": "3.10.2"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "a42ccb73e7d9bfdf27e036f1d2b8b681e55fc0743cc5586bc2474d4a60f4b886"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
